Hands-on Exercise 3: 1st and 2nd Order Spatial Point Patterns Analysis (SPPA) Methods

Author

William

Published

August 30, 2024

Modified

August 31, 2024

Overview

Lesson Notes

Distribution of Crime Incidents
  • Point pattern analysis (PPA) is the study of the spatial arrangements of points in (usually 2D) space to find out if the pattern is random or the result of something.

  • Spatial PPA Methods

    • 1st order

      • Density-based

        • Kernel density estimation: compute intensity of a point distribution (Adaptive v. Fixed Bandwidth)

        • Quadrat analysis to obtain variance-mean ratio (VMR: >0,0,<0): measure of dispersion rather than a measure of pattern; single measure for the entire distribution, so variation within the region are not recognised

      • Distance-based

        • Nearest Neighbour Index

        • G-function

        • F-function

        • K-function

        • L-function

Introduction

Using appropriate functions of spatstat, this hands-on exercise aims to discover the spatial point processes of childecare centres in Singapore.

The specific questions we would like to answer are as follows:

  • are the childcare centres in Singapore randomly distributed throughout the country?

  • if the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?

Data

  • CHILDCARE, a point feature data providing both location and attribute information of childcare centres. It was downloaded from Data.gov.sg and is in geojson format.

  • MP14_SUBZONE_WEB_PL, a polygon feature data providing information of URA 2014 Master Plan Planning Subzone boundary data. It is in ESRI shapefile format. This data set was also downloaded from Data.gov.sg.

  • CostalOutline, a polygon feature data showing the national boundary of Singapore. It is provided by SLA and is in ESRI shapefile format.

Packages

  • sf, a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.

  • spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.

  • raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.

  • maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.

  • tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)

Spatial Data Wrangling

Import the 3 geospatial datasets.

childcare_sf <- st_read('data/child-care-services-geojson.geojson') %>% st_transform(crs = 3414)
Reading layer `child-care-services-geojson' from data source 
  `/Users/williamtjw/is415-gaa-williamtjw/resources/Hands-on_Ex/ex3/data/child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

Note that the CRS for childcare_sf is WGS84, not SVY21.

childcare_sf3414 <- st_transform(childcare_sf, crs = 3414)

To check,

st_geometry(childcare_sf3414)
Geometry set for 1545 features 
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 5 geometries:
POINT Z (27976.73 45716.7 0)
POINT Z (25824 29900.09 0)
POINT Z (31399.04 37416.36 0)
POINT Z (29268.43 40942.1 0)
POINT Z (41217.74 33554.94 0)
sg_sf <- st_read(dsn = 'data/', layer = 'CostalOutline')
Reading layer `CostalOutline' from data source 
  `/Users/williamtjw/is415-gaa-williamtjw/resources/Hands-on_Ex/ex3/data' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21
mpsz_sf <- st_read(dsn = 'data/', layer = 'MP14_SUBZONE_WEB_PL')
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `/Users/williamtjw/is415-gaa-williamtjw/resources/Hands-on_Ex/ex3/data' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Ensure that they are projected in same projection system

st_crs(childcare_sf3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Note that even though mpsz_sf and sg_sf are projected in svy21, the EPSG codes are 9001 instead of 3414.

Use st_set_crs() of sf package to assign the correct EPSG code

sg_sf3414 <- st_set_crs(sg_sf, 3414)
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that
mpsz_sf3414 <- st_set_crs(mpsz_sf, 3414)
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that

To check,

st_crs(sg_sf3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(mpsz_sf3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_geometry(sg_sf3414)
Geometry set for 60 features 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21 / Singapore TM
First 5 geometries:
POLYGON ((14362.86 32307.49, 14414.06 32336.75,...
POLYGON ((25683.97 26236.91, 25696.79 26245.47,...
POLYGON ((11471.97 46273.01, 11483.97 46279.01,...
POLYGON ((12602.3 32061.35, 12622.55 32066.42, ...
POLYGON ((17915.53 46770.73, 17933.13 46772.76,...
st_geometry(mpsz_sf3414)
Geometry set for 323 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM
First 5 geometries:
MULTIPOLYGON (((31495.56 30140.01, 31980.96 296...
MULTIPOLYGON (((29092.28 30021.89, 29119.64 300...
MULTIPOLYGON (((29932.33 29879.12, 29947.32 298...
MULTIPOLYGON (((27131.28 30059.73, 27088.33 297...
MULTIPOLYGON (((26451.03 30396.46, 26440.47 303...

Plotting to show spatial patterns

tmap_mode('plot')
tmap mode set to plotting
tm_shape(mpsz_sf3414) + tm_polygons() +
  tm_shape(childcare_sf3414) + tm_dots()

Notice that all the geospatial layers are within the same map extend. This shows that their referencing system and coordinate values are referred to similar spatial context. This is very important in any geospatial analysis.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf)+
  tm_dots()

This is the default background of the internet map layer: ESRI.WorldGrayCanvas. Currently, three internet map layers: ESRI.WorldGrayCanvas, OpenStreetMap, and ESRI.WorldTopoMap.

tmap_mode('plot')
tmap mode set to plotting

Always remember to switch back to plot mode after the interactive map. This is because, each interactive mode will consume a connection. You should also avoid displaying ecessive numbers of interactive maps (i.e. not more than 10) in one RMarkdown document when publish on Netlify.

Geospatial Data wrangling

Many geospatial analysis packages require the input geospatial data in sp’s Spatial* classes.

Converting sf data frames to sp’s Spatial* class

childcare <- as_Spatial(childcare_sf3414)
sg <- as_Spatial(sg_sf3414)
mpsz <- as_Spatial(mpsz_sf3414)

Examine the information of these three Spatial* classes.

structure(childcare) # class: SpatialPointsDataFrame
class       : SpatialPointsDataFrame 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                  <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
structure(sg) # class: SpatialPolygonsDataFrame
class       : SpatialPolygonsDataFrame 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       : GDO_GID, MSLINK, MAPID,              COSTAL_NAM 
min values  :       1,      1,     0,             ISLAND LINK 
max values  :      60,     67,     0, SINGAPORE - MAIN ISLAND 
structure(mpsz) # class: SpatialPolygonsDataFrame
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 

Converting the Spatial* class into generic sp format

childcare_sp <- as(childcare, 'SpatialPoints')
sg_sp <- as(sg, 'SpatialPolygons')
mpsz_sp <- as(mpsz, 'SpatialPolygons')
childcare_sp
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sg_sp
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
mpsz_sp
class       : SpatialPolygons 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Converting the generic sp format into spatstat’s ppp format

Use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format

childcare_ppp <- as.ppp(childcare_sf)
Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
marks
childcare_ppp
Marked planar point pattern: 1545 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

Plot childcare_ppp

plot(childcare_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

Examine the ppp object’s summary stats.

summary(childcare_ppp)
Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

Notice the warning message about duplicates. In spatial point patterns analysis an issue of significant is the presence of duplicates. The statistical methodology used for spatial point patterns processes is based largely on the assumption that process are simple, that is, that the points cannot be coincident.

Handling duplicated points

Check for duplicates in the ppp object.

any(duplicated(childcare_ppp))
[1] FALSE

To count the number of co-incidience points,

# multiplicity(childcare_ppp)

To count how many locations have more than one point event,

sum(multiplicity(childcare_ppp) > 1)
[1] 0

View the locations of these duplicate point events using childcare.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) +
  tm_dots(alpha=0.4, 
          size=0.05)

Darker spots indicate duplicate point events

tmap_mode('plot')
tmap mode set to plotting

Solutions

  • The easiest way is to delete the duplicates. But, that will also mean that some useful point events will be lost.

  • The second solution is use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

  • The third solution is to make each point “unique” and then attach the duplicates of the points to the patterns as marks, as attributes of the points. Then you would need analytical techniques that take into account these marks.

Jittering approach:

childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)
# rjitter applies independent random displacements to each point in a point pattern.
sum(multiplicity(childcare_ppp_jit)>1)
[1] 0

Creating owin object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.n

Covert sg SpatialPolygon object into owin object of spatstat.

sg_owin <- as.owin(sg_sf)
plot(sg_owin)

summary(sg_owin)
Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

Combining point events object and owin object

Extract childcare events that are located within Singapore

childcareSG_ppp = childcare_ppp[sg_owin]
summary(childcareSG_ppp)
Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401
plot(childcareSG_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

1st order SPPA using spatstat

  • deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes,

  • performing Confirmatory SPPA by using Nearest Neighbour (NN) statistics.

KDE (of childcare services in Singapore)

Using automatic bandwidth selection method

Compute a kde by using the following configurations of density() of spatstat:

  • bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl().

  • The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”.

  • The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.

kde_childcareSG_bw <- density(
  childcareSG_ppp,
  sigma=bw.diggle,
  edge=TRUE,
  kernel='gaussian')

# 'density()' computes kernel density estimates. Its default method does so with the given kernel and bandwidth for univariate observations.

# 'sigma' computes the covariance matrix of errors (bandwidth of kernel)

# 'edge = TRUE': This parameter adjusts the density estimation near the boundaries to reduce edge effects, which can cause underestimation of densities near the boundaries of the study area
plot(kde_childcareSG_bw)

The density values of the output range from 0 to 0.000035. This is because the default unit of measurement of svy21 is in meter. As a result, the density values computed is in “number of points per square meter”.

Retrieve the bandwidth (sigma) used to compute the kde layer.

bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
298.4095 

Rescalling KDE values

Use rescale.ppp() to convert the unit of measurement from m to km

childcareSG_ppp_km <- rescale.ppp(childcareSG_ppp,1000,'km')

Re-run density() using the rescaled data and plot the kde map

kde_childcareSG_bw <- density(
  childcareSG_ppp_km,
  sigma=bw.diggle,
  edge=TRUE,
  kernel='gaussian'
)
plot(kde_childcareSG_bw)

KDE map is similar to the first but with a different scale in the legend

Working with different automatic bandwidth methods

Beside bw.diggle(), there are three other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().

bw.CvL(childcareSG_ppp_km)
   sigma 
4.543278 
bw.scott(childcareSG_ppp_km)
 sigma.x  sigma.y 
2.224898 1.450966 
bw.ppl(childcareSG_ppp_km)
    sigma 
0.3897114 
bw.diggle(childcareSG_ppp_km)
    sigma 
0.2984095 

Baddeley et. (2016):

  • bw.ppl() tends to produce more appropriate values when the pattern consists predominantly of tight clusters

  • if the purpose of once study is to detect a single tight cluster in the midst of random noise then the bw.diggle() method seems to work best

bw.diggle() v. bw.ppl() output

kde_childcareSG_ppl <- density(
  childcareSG_ppp_km,
  sigma=bw.ppl,
  edge=TRUE,
  kernel='gaussian'
)

par(mfrow=c(1,2))
plot(kde_childcareSG_bw,main='bw.diggle')
plot(kde_childcareSG_ppl,main='bw.ppl')

Working with different kernel methods

By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics.

par(mfrow=c(2,2))
plot(density(childcareSG_ppp_km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp_km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
Warning in density.ppp(childcareSG_ppp_km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp_km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
Warning in density.ppp(childcareSG_ppp_km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp_km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")
Warning in density.ppp(childcareSG_ppp_km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

# If your editor window is too small, might get the error: "Error in plot.new() : figure margins too large".
# By writing 
par(mar=c(1,1,1,1)) 
#instead, the error is fixed; or u could simply resize your window.

Fixed and Adaptive KDE

Compute KDE using fixed BW

kde_childcareSG_600 <- density(
  childcareSG_ppp_km,
  sigma=0.6,
  edge=TRUE,
  kernel='gaussian'
)

# sigma = 0.6 instead of 600 because
# childcareSG_ppp_km object is in KILOMETERS
plot(kde_childcareSG_600)

Compute using adaptive BW

Fixed bandwidth method is very sensitive to highly skewed distribution of spatial point patterns over geographical units (for example urban versus rural). One solution is to use adaptive BW instead.

Derive adaptive KDE using density.adaptive() of spatstat.

kde_childcareSG_adaptive <- adaptive.density(
  childcareSG_ppp_km,
  method = 'kernel' # voronoi, kernel, nearest being the other methods
)
plot(kde_childcareSG_adaptive)

Compare the fixed and adaptive KDE outputs.

par(mfrow=c(1,2))
plot(kde_childcareSG_bw, main = "Fixed BW")
plot(kde_childcareSG_adaptive, main = "Adaptive BW")

Converting KDE output into grid object

Convert it so that it is suitable for mapping purposes.

gridded_kde_childcareSG_bw <- as(kde_childcareSG_bw,'SpatialGridDataFrame')
spplot(gridded_kde_childcareSG_bw)

Converting gridded output into raster

Convert the gridded kernal density objects into RasterLayer object by using raster() of raster package

kde_childcareSG_bw_raster <- raster(kde_childcareSG_bw)
kde_childcareSG_bw_raster # min value differs from reference output(?)
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -1.005814e-14, 28.51831  (min, max)

Assigning projection systems

Include the CRS information on kde_childcareSG_bw_raster RasterLayer

projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : -1.005814e-14, 28.51831  (min, max)

Note the changes in the ‘CRS’ property.

Visualising the output in tmap

Display the raster in cartographic quality map

tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("layer", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

Comparing Spatial Point Patterns using KDE

Compare KDE of childcare at Punggol, Tampines, Chua Chu Kang and Jurong West planning areas.

Extract and Plot study areas

pg <- mpsz_sf %>% filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>% filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>% filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>% filter(PLN_AREA_N == "JURONG WEST")
par(mfrow=c(2,2))
plot(pg, main = "Punggol")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(tm, main = "Tampines")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(ck, main = "Choa Chu Kang")
Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

plot(jw, main = "Jurong West")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

Create owin object

  • convert these sf objects into owin objects that is required by spatstat.
pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)

Combine childcare points and the study area

  • extract childcare that is within the specific region for future analysis
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]
  • convert m to km using rescale.ppp()
childcare_pg_ppp_km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp_km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp_km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp_km = rescale.ppp(childcare_jw_ppp, 1000, "km")
  • plot these four study areas and the locations of the childcare centres
par(mfrow=c(2,2))
plot(childcare_pg_ppp_km, main="Punggol")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_tm_ppp_km, main="Tampines")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 89 symbols are shown in the symbol map
plot(childcare_ck_ppp_km, main="Choa Chu Kang")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_jw_ppp_km, main="Jurong West")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 88 symbols are shown in the symbol map

  • compute KDE (bw.diggle method is used to derive each of the bandwidth)
par(mfrow=c(2,2))
plot(density(childcare_pg_ppp_km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp_km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")
plot(density(childcare_ck_ppp_km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp_km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")

  • Compute fixed bandwidth KDE (250m for comparison)
par(mfrow=c(2,2))
plot(density(childcare_ck_ppp_km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp_km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")
plot(density(childcare_pg_ppp_km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp_km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

Nearest Neighbour Analysis

Perform the Clark-Evans test of aggregation for a spatial point pattern by using clarkevans.test() of statspat.

  • Ho = The distribution of childcare services are randomly distributed.

  • H1= The distribution of childcare services are not randomly distributed.

  • 95% confidence level

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"), # type of alternative for the hypothesis test
                nsim=99                     # Number of Monte Carlo simulations to perform
  )

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
# clarkevans test is a crude measure of clustering or ordering of a point pattern.

Strong statistical support to accept H1. It makes sense as much resources and planning are put into such publics services so their locations are very intentional and unlikely to be randomly built across Singapore.

Clark and Evans Test: Choa Chu Kang planning area

clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.90613, p-value = 0.1608
alternative hypothesis: two-sided

Clark and Evans Test: Tampines planning area

clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.75977, p-value = 1.453e-05
alternative hypothesis: two-sided

References